Bulk analysis of MS lesions(Supervised)

This is an analysis performed on the bulk RNA-Seq dataset

The workflow is as follows:
1. Inspect variance in the dataset (Filtering, Normalization, PCA)
2. Perform differential expression analysis (Using limma)
3. Annotate biological signatures found by DEA (GSEA)

library(edgeR)
library(fgsea)
library(qusage)
library(DESeq2)
library(org.Hs.eg.db)
library(ggplot2)
library(enrichR)
library(VennDiagram)
library(ggpubr)
library(ggrepel)
library(PCAtools)
library(RColorBrewer)
library(ggiraph)
library(UpSetR)
library(rhandsontable)
library(nloptr)
library(ggVennDiagram)
library(BiocManager)
library(stringr)
library(knitr)

Counts list

dat = read.csv('/Users/skander/Library/CloudStorage/OneDrive-Personnel/Documents/UDEM/Master/BulkData/data/bulkExtra/Totalbulk/merged_counts.csv',sep=';', header = TRUE, row.names=1)
knitr::kable(dat[1:12, 1:12], caption ="A Knitr kable")
gene_name S13_b1 S14_b1 S23 S26 S29 S33 S37 S38 S3_b1 S43 S45
ENSG00000000003 TSPAN6 131.00 46.00 83.00 59.00 79.00 98.00 56.00 57.00 220.00 74.00 725.00
ENSG00000000005 TNMD 1.00 1.00 0.00 2.00 2.00 1.00 0.00 0.00 5.00 1.00 3.00
ENSG00000000419 DPM1 81.00 103.00 121.00 128.00 141.00 91.00 61.00 68.00 66.00 122.00 112.00
ENSG00000000457 SCYL3 91.24 105.74 120.97 107.41 115.13 93.44 64.38 69.50 74.44 168.65 191.43
ENSG00000000460 C1orf112 95.86 56.26 72.03 68.59 67.87 44.56 59.62 50.50 91.56 88.35 163.57
ENSG00000000938 FGR 28.00 11.00 15.00 16.00 6.00 58.00 44.00 28.00 56.00 16.00 113.00
ENSG00000000971 CFH 103.00 41.00 90.00 70.00 41.00 198.00 59.00 151.00 100.00 97.00 105.00
ENSG00000001036 FUCA2 104.60 70.29 65.07 75.40 73.47 83.00 62.29 71.67 89.90 81.02 228.64
ENSG00000001084 GCLC 482.00 451.00 406.00 442.00 500.00 466.97 395.00 344.00 388.00 559.00 1042.00
ENSG00000001167 NFYA 533.00 325.00 224.00 333.97 307.00 340.00 508.89 344.00 450.00 260.00 491.00
ENSG00000001460 STPG1 173.02 227.16 157.14 155.07 183.09 140.37 131.45 146.10 182.06 205.00 214.00
ENSG00000001461 NIPAL3 809.98 963.84 656.86 670.93 757.91 776.63 1438.55 675.90 692.94 1072.00 992.00

A Knitr kable

Meta data

meta = read.csv('/Users/skander/Library/CloudStorage/OneDrive-Personnel/Documents/UDEM/Master/BulkData/data/bulkExtra/Totalbulk/meta_data_3batch.csv', header = T,sep = ';', row.names = 1)
names(meta)[names(meta) == "type"] <- "lesion_type"
meta$sex = 'M'
meta[which(meta$sample=="AB203"),]$sex = 'F'
meta[which(meta$sample=="AB253"),]$sex = 'F'

meta$disease="MS"
meta[which(meta$sample %in% c("AB257","AB251","AB252")),]$disease="HC"
#meta$region="WM"
meta$group=paste(meta$disease,meta$tissue,sep=".")
meta$reads=colSums(dat[-1])[rownames(meta)]

knitr::kable(meta, caption ="A Knitr kable")
id batch lesion_type tissue condition sample typeL sex disease group reads
S13_b1 S13_b1 batch1 WML WM Lesion AB187 M MS MS.WM 10675040
S14_b1 S14_b1 batch1 NAGM Cortical Normal AB187 M MS MS.Cortical 11601256
S23 S23 batch1 CL Cortical Lesion AB203 F MS MS.Cortical 10388672
S26 S26 batch1 WML WM Lesion AB187 M MS MS.WM 11312959
S29 S29 batch1 CL Cortical Lesion AB187 M MS MS.Cortical 10849903
S33 S33 batch1 CL Cortical Lesion AB200 M MS MS.Cortical 15027141
S37 S37 batch1 WML WM Lesion AB200 M MS MS.WM 12031415
S38 S38 batch1 NAGM Cortical Normal AB200 M MS MS.Cortical 12856284
S3_b1 S3_b1 batch1 PVML PV Lesion AB200 M MS MS.PV 10431353
S43 S43 batch1 CL Cortical Lesion AB203 F MS MS.Cortical 12121099
S45 S45 batch1 PVML PV Lesion AB203 F MS MS.PV 10093676
S53 S53 batch1 NAWM WM Normal AB187 M MS MS.WM 10907244
S55 S55 batch1 NAGM Cortical Normal AB200 M MS MS.Cortical 20126720
S62 S62 batch1 WML WM Lesion AB200 M MS MS.WM 12859996
S64 S64 batch1 NAGM Cortical Normal AB203 F MS MS.Cortical 13502425
S67 S67 batch1 WML WM Lesion AB203 F MS MS.WM 11484555
S70 S70 batch1 WML WM Lesion AB203 F MS MS.WM 10809418
S71 S71 batch1 NAWM WM Normal AB187 M MS MS.WM 13315868
S72 S72 batch1 CL Cortical Lesion AB187 M MS MS.Cortical 12938331
S76 S76 batch1 CL Cortical Lesion AB200 M MS MS.Cortical 14486529
S77 S77 batch1 WML WM Lesion AB200 M MS MS.WM 14583455
S79 S79 batch1 NAWM WM Normal AB203 F MS MS.WM 14931094
S9_b1 S9_b1 batch1 NAGM Cortical Normal AB187 M MS MS.Cortical 14719787
Sample1 Sample1 batch2 PVML PV Lesion AB200 M MS MS.PV 20735000
Sample10 Sample10 batch2 CL Cortical Lesion AB187 M MS MS.Cortical 9925960
Sample11 Sample11 batch2 WML WM Lesion AB187 M MS MS.WM 19113957
Sample12 Sample12 batch2 CL Cortical Lesion AB187 M MS MS.Cortical 13267587
Sample15 Sample15 batch2 WML WM Lesion AB187 M MS MS.WM 14887133
Sample16 Sample16 batch2 WML WM Lesion AB187 M MS MS.WM 12598401
Sample17 Sample17 batch2 PVML PV Lesion AB203 F MS MS.PV 10162274
Sample18 Sample18 batch2 CL Cortical Lesion AB203 F MS MS.Cortical 11609221
Sample21 Sample21 batch2 NAWM WM Normal AB203 F MS MS.WM 11016594
Sample24 Sample24 batch2 NAWM WM Normal AB203 F MS MS.WM 14238903
Sample25 Sample25 batch2 NAGM Cortical Normal AB187 M MS MS.Cortical 15168403
Sample34 Sample34 batch2 WML WM Lesion AB200 M MS MS.WM 18291181
Sample46 Sample46 batch2 NAGM Cortical Normal AB203 F MS MS.Cortical 16625491
Sample52 Sample52 batch2 NAWM WM Normal AB187 M MS MS.WM 12701779
Sample56 Sample56 batch2 PVML PV Lesion AB200 M MS MS.PV 13165303
Sample57 Sample57 batch2 CL Cortical Lesion AB200 M MS MS.Cortical 11305279
Sample6 Sample6 batch2 NAGM Cortical Normal AB200 M MS MS.Cortical 22662880
Sample60 Sample60 batch2 NAWM WM Normal AB200 M MS MS.WM 14547748
Sample65 Sample65 batch2 PVML PV Lesion AB203 F MS MS.PV 13327909
Sample66 Sample66 batch2 WML WM Lesion AB203 F MS MS.WM 8626915
Sample7 Sample7 batch2 WML WM Lesion AB200 M MS MS.WM 17559199
Sample74 Sample74 batch2 NAGM Cortical Normal AB187 M MS MS.Cortical 17633780
Sample75 Sample75 batch2 NAWM WM Normal AB200 M MS MS.WM 20237852
S1 S1 batch3 GM Cortical AB257 GM M HC HC.Cortical 14746041
S2 S2 batch3 GM Cortical AB257 GM M HC HC.Cortical 16135512
S3_b3 S3_b3 batch3 WM WM AB257 WM M HC HC.WM 17280964
S4 S4 batch3 WM WM AB257 WM M HC HC.WM 16128150
S5 S5 batch3 PVM PV AB257 PVM M HC HC.PV 15969401
S6 S6 batch3 PVM PV AB257 PVM M HC HC.PV 14119977
S7 S7 batch3 GM Cortical AB252 GM M HC HC.Cortical 16993959
S8 S8 batch3 GM Cortical AB252 GM M HC HC.Cortical 16902227
S9_b3 S9_b3 batch3 WM WM AB252 WM M HC HC.WM 14200304
S10 S10 batch3 WM WM AB252 WM M HC HC.WM 14993001
S11 S11 batch3 PVM PV AB252 PVM M HC HC.PV 13322475
S12 S12 batch3 PVM PV AB252 PVM M HC HC.PV 11873861
S13_b3 S13_b3 batch3 GM Cortical AB251 GM M HC HC.Cortical 14540503
S14_b3 S14_b3 batch3 GM Cortical AB251 GM M HC HC.Cortical 15068466
S15 S15 batch3 WM WM AB251 WM M HC HC.WM 13301374
S16 S16 batch3 WM WM AB251 WM M HC HC.WM 15610086
S17 S17 batch3 PVM PV AB251 PVM M HC HC.PV 12503687
S18 S18 batch3 PVM PV AB251 PVM M HC HC.PV 15465598
S102f S102f batch3 CL Cortical AB253 CL F MS MS.Cortical 11101420
S104f S104f batch3 CL Cortical AB253 CL F MS MS.Cortical 16624705
S106f S106f batch3 CL Cortical AB253 CL F MS MS.Cortical 18467343
S108f S108f batch3 NAGM Cortical AB253 NAGM F MS MS.Cortical 12319237
S109f S109f batch3 NAGM Cortical AB253 NAGM F MS MS.Cortical 14185427
S110f S110f batch3 CL Cortical AB253 CL F MS MS.Cortical 12999955
S112f S112f batch3 NAWM WM AB253 NAWM F MS MS.WM 13363398
S113f S113f batch3 NAWM WM AB253 NAWM F MS MS.WM 16157577
S114f S114f batch3 CL Cortical AB253 CL F MS MS.Cortical 16266665
S115f S115f batch3 NAGM Cortical AB253 NAGM F MS MS.Cortical 17885078
S118f S118f batch3 PVML PV AB253 Mixed.PD F MS MS.PV 10443427
S119f S119f batch3 PVML PV AB253 Active.PD F MS MS.PV 13248535
S120f S120f batch3 PVML PV AB253 Active.PD F MS MS.PV 12139042
S121f S121f batch3 PVML PV AB253 Mixed.PD F MS MS.PV 10019098
S19a S19a batch3 WML WM AB253 Active F MS MS.WM 17274976
S22a S22a batch3 PVML PV AB253 PVML F MS MS.PV 11231127
S28b S28b batch3 WML WM AB187 Active.PD M MS MS.WM 15599429
S30b S30b batch3 NAWM WM AB187 NAWM M MS MS.WM 16356084
S31b S31b batch3 WML WM AB187 Active.PD M MS MS.WM 13453907
S51c S51c batch3 CL Cortical AB187 CL M MS MS.Cortical 16757776
S69c S69c batch3 PVML PV AB203 PVML F MS MS.PV 13814668
S73d S73d batch3 WML WM AB187 Active.PD M MS MS.WM 13243821
S78d S78d batch3 NAWM WM AB200 NAWM M MS MS.WM 14215553
S81d S81d batch3 PVML PV AB203 PVML F MS MS.PV 15303502
S89f S89f batch3 WML WM AB253 Early.Active F MS MS.WM 10392586
S90f S90f batch3 WML WM AB253 Early.Active F MS MS.WM 18177867
S91f S91f batch3 WML WM AB253 Early.Active F MS MS.WM 11856582
S92f S92f batch3 NAGM Cortical AB253 NAGM F MS MS.Cortical 12412586
S94f S94f batch3 WML WM AB253 Mixed.PD F MS MS.WM 13343462
S96f S96f batch3 WML WM AB253 Active.AD F MS MS.WM 8964589

A Knitr kable

Define functions for the analysis: - find_enrichments(): For a given comparison (DEA) extract enrichments for GO(BP), GO(MF), KEGG and TFT
- fisher_enrichment(): Perform an enrichment of genesets using a fisher test approach - plot.gene(): Plot gene expression for a given gene across conditions.

The gmt files contain the genesets for GSEA

go = read.gmt('/Users/skander/Library/CloudStorage/OneDrive-Personnel/Documents/UDEM/Master/BulkData/GO/c5.go.bp.v7.5.1.symbols.gmt')
mf = read.gmt('/Users/skander/Library/CloudStorage/OneDrive-Personnel/Documents/UDEM/Master/BulkData/GO/c5.go.mf.v7.5.1.symbols.gmt')
kegg = read.gmt('/Users/skander/Library/CloudStorage/OneDrive-Personnel/Documents/UDEM/Master/BulkData/GO/c2.cp.kegg.v7.5.1.symbols.gmt')
tft = read.gmt('/Users/skander/Library/CloudStorage/OneDrive-Personnel/Documents/UDEM/Master/BulkData/GO/c3.tft.v7.5.1.symbols.gmt')

sign_cols=c(DN='deepskyblue',NS='grey80',UP='tomato')

plot.gene <- function(x){
  gene_id = rownames(subset(resA,gene_name==x))[1]
  meta$gene = y$E[gene_id,rownames(meta)]
  ggplot(meta,aes(lesion_type,gene))+geom_boxplot()+geom_point(aes(col=batch))+theme_bw()+ ggtitle(x)
}

fisher_enrichment<-function(cluster_markers,pathway,universe,pathway_name){
enrichment_tables=list() 
  for(cluster in 1:length(cluster_markers)){
    cluster_name=names(cluster_markers)[cluster]
    output <- lapply(pathway, function(x) {
    freq.table <- table(factor(universe %in% as.character(cluster_markers[[cluster]]), 
                          levels = c(TRUE,FALSE)), 
                          factor(universe %in% x, 
                          levels = c(TRUE,FALSE)))

      fit <- fisher.test(freq.table, alternative = "greater")
      interSection <- intersect(cluster_markers[[cluster]], x)
      interSection <- paste(interSection, collapse = ",")
      return(value = c(NOM_pval = fit$p.value, INTERSECT = interSection,"OR"=fit$estimate))})
    
    term_names=character()
    for (pathway.term in 1:length(output)){
      term_names[pathway.term]=pathway[[pathway.term]][1]
    }
  
    results=data.frame(do.call("rbind",output))
    results$fdr=p.adjust(as.numeric(as.character(results$NOM_pval)),method = "BH")
    results=results[order(results$fdr),]
    enrichment_tables[[cluster_name]]=results
  }

  return(enrichment_tables)
}

find_enrichments<-function(resA,n=5,show='both'){
  resA.t = resA$t
  names(resA.t)=resA$gene_name 
  if(show=='up'){
    resA.bp = subset(fgsea(go,resA.t),NES>0)
    resA.mf = subset(fgsea(mf,resA.t),NES>0)
    resA.kegg= subset(fgsea(kegg,resA.t),NES>0)
    resA.tft= subset(fgsea(tft,resA.t),NES>0)
  }else if(show=='down'){
    resA.bp = subset(fgsea(go,resA.t),NES<0)
    resA.mf = subset(fgsea(mf,resA.t),NES<0)
    resA.kegg= subset(fgsea(kegg,resA.t),NES<0)   
    resA.tft= subset(fgsea(tft,resA.t),NES<0)
    
  }else{
    resA.bp = fgsea(go,resA.t)
    resA.mf = fgsea(mf,resA.t)
    resA.kegg= fgsea(kegg,resA.t)   
    resA.tft= subset(fgsea(tft,resA.t))
  }
  
  resA.bp$GO = 'BP'
  resA.mf$GO = 'MF'
  resA.kegg$GO = 'KEGG'
  out=list(enr=list(bp=resA.bp,mf=resA.mf,kegg=resA.kegg,tft=resA.tft))
  
  enr = rbind(
    head(resA.bp[order(resA.bp$padj),],n),
    head(resA.mf[order(resA.mf$padj),],n),
    head(resA.kegg[order(resA.kegg$padj),],n)
  )
  enr=data.frame(enr)
  enr$pathway=gsub('_',' ',gsub('GO_|KEGG_','',enr$pathway))
  enr$pathway=factor(enr$pathway,levels=as.character(enr[order(enr$GO,enr$NES),'pathway']))
  p=ggplot(enr,aes(pathway,NES,fill=GO))+geom_col(width=.6,alpha=.6)+geom_col(width=.6,fill=NA,aes(col=GO))+coord_flip()+theme_bw()+geom_vline(xintercept=c(n,n*2)+.5,linetype='dashed')+geom_hline(yintercept=0)
  out$plot=p
  return(out)
}

##top pathway
topPath<- function(fish ,specific, p ,n){
  l=n*20
  result <- matrix(NA, l, l)
  
  colnames(result) <- rownames(result) <- rownames(getElement(fish, specific))[1:l]
  for(i in 1:l){
    for(j in 1:l){
      minlist = 0
      if(length(unlist(strsplit(fish[[specific]][["INTERSECT"]][i], split = ",")))>=length(unlist(strsplit(fish[[specific]][["INTERSECT"]][j], split = ",")))) minlist=length(unlist(strsplit(fish[[specific]][["INTERSECT"]][j], split = ","))) 
      else minlist=length(unlist(strsplit(fish[[specific]][["INTERSECT"]][i], split = ",")))
      result[i, j] <- length(intersect(unlist(strsplit(fish[[specific]][["INTERSECT"]][i], split = ",")), unlist(strsplit(fish[[specific]][["INTERSECT"]][j], split = ","))))/(minlist)      
      result[j, i] <- result[i, j]     
    }
  }
  listBPspecific = list(c(rownames(result)[1]))
  i=1
  while(i < l){
    j=i
    while(j < l){
      if(result[i, j] < p){ 
        i=j
        listBPspecific <- c(listBPspecific, rownames(result)[j])
        j=l
      }else {
        j=j+1
      }
    }
    if(length(listBPspecific)>n-1) 
      i=l
  }
  return(listBPspecific)
}

Counts list et metadata for batch3

datB3 <- dat[,48:95]
metaB3 <- meta[47:94,]

Differential expression analysis with limma-voom

Filtering to remove lowly expressed genes + Normalization for composition bias

use the edgeR/limma workflow to filter out lowly expressed genes. Voom is used to normalize the raw counts. A linear mixed model is used to fit the expression values. Defining the model is done with the model.matrix function of edgeR. In the model are included factors of interest (group). If a batch/patient effect is detected, it can be added to the model to account for its effect.

Only Batch3

Without filtering low expressed genes

d0 <- DGEList(datB3)
d0 <- calcNormFactors(d0)
cutoff <- 5
drop <- which(apply(cpm(d0), 1, max) < cutoff)
d <- d0[-drop,]

mmB3 <- model.matrix(~ 0 + lesion_type+sex+sample+reads, data=metaB3[colnames(d),])

y <- voom(d0[,rownames(metaB3)], mmB3, plot = T)
## Coefficients not estimable: sampleAB253 sampleAB257

Filtering low expressed genes

# Filtering low expressed genes
yB3 <- voom(d[,rownames(metaB3)], mmB3, plot = T)
## Coefficients not estimable: sampleAB253 sampleAB257

ALL batchs

Without filtering low expressed genes

d0 <- DGEList(dat[-1])
d0 <- calcNormFactors(d0)
cutoff <- 5
drop <- which(apply(cpm(d0), 1, max) < cutoff)
d <- d0[-drop,]

mm <- model.matrix(~ 0 + lesion_type+batch+sample+sex, data=meta[colnames(d),])

y <- voom(d0[,rownames(meta)], mm, plot = T)
## Coefficients not estimable: sampleAB257 sexM

Filtering low expressed genes

# Filtering low expressed genes
y <- voom(d[,rownames(meta)], mm, plot = T)
## Coefficients not estimable: sampleAB257 sexM

Reads By Batch

meta$id = factor(meta$id, levels=meta[order(meta$sample),'id'])
ggplot(meta,aes(id, reads, fill=sample))+theme_bw()+geom_col(width=0.6)+ theme(axis.text.x=element_text(angle=60,hjust=1,size=5))+facet_grid(~batch, space = "free_x",scales = "free_x")

PCA Batch3

Perform a PCA on the normalized samples and identify main factors of variation.

col by Sample

pcaB3<-pca(yB3$E, metadata = metaB3)
pcaB3S <- biplot(pcaB3, colby = 'sample',hline = 0, vline = 0,legendPosition = 'right', lab="",title = "")
pcaB3S

### col by lesion type

pcaB3<-pca(yB3$E, metadata = metaB3)
pcaB3L <- biplot(pcaB3, colby = 'lesion_type',hline = 0, vline = 0,legendPosition = 'right', lab="",title = "")
pcaB3L

### col by tissue

pcaB3<-pca(yB3$E, metadata = metaB3)
pcaB3t <- biplot(pcaB3, colby = 'tissue',hline = 0, vline = 0,legendPosition = 'right', lab="",title = "")
pcaB3t

### col by type

pcaB3<-pca(yB3$E, metadata = metaB3)
pcaB3tl <- biplot(pcaB3, colby = 'typeL',shape='tissue',hline = 0, vline = 0,legendPosition = 'right', lab="",title = "")
pcaB3tl

PCA All DATA

Perform a PCA on the normalized samples and identify main factors of variation.

Data including batch effect

BD<-pca(y$E, metadata = meta)
BD <- biplot(BD, colby = 'batch',hline = 0, vline = 0,legendPosition = 'right', lab="",title = "Data including batch effect")
BD

BDCorr <- removeBatchEffect(y$E,batch = meta$batch)

Batch effect removed from data

BDCorr <- removeBatchEffect(y$E,batch = meta$batch)
BD1 <- pca(BDCorr, metadata = meta)
BDb <- biplot(BD1, colby = 'batch',hline = 0, vline = 0,legendPosition = 'right', lab="",title="Batch effect removed from data")
BDb

Effect of samples on data (sex, type of ms,…) (col by Sample)

BDsample <- biplot(BD1, colby = 'sample',hline = 0, vline = 0,legendPosition = 'right', lab="",title="Effect of samples on data (sex, type of ms,...)")
BDsample

Effect of samples on data (sex, type of ms,…) (col by Sex)

BDsex <- biplot(BD1, colby = 'sex',hline = 0, vline = 0,legendPosition = 'right', lab="",title="Effect of samples on data (sex, type of ms,...)")
BDsex

Batch+Patient effect removed from data (col by Sex)

BDCorrBS <- removeBatchEffect(y$E,batch = meta$batch, batch2 = meta$sample)
BD2 <- pca(BDCorrBS, metadata = meta)
BDsexCorr <- biplot(BD2, colby = 'sex',hline = 0, vline = 0,legendPosition = 'right', lab="",title="Batch+Patient effect removed from data")
BDsexCorr

Batch+Patient effect removed from data (col by Sample)

BDsampleCorr <- biplot(BD2, colby = 'sample',hline = 0, vline = 0,legendPosition = 'right', lab="",title="Batch+Patient effect removed from data")
BDsampleCorr

Batch+Patient effect removed from data (col by Lesion type and Sex)

BDg <- biplot(BD2, colby = 'lesion_type',shape = 'sex',hline = 0, vline = 0,legendPosition = 'right', lab="",title="PCA")
## PCA
BDg

gene name reference(ref_G)

conv_B=AnnotationDbi::select(org.Hs.eg.db,keys=rownames(y$E),keytype='ENSEMBL',columns=c('GENENAME','SYMBOL'))
ref_G=conv_B[!duplicated(conv_B[,1]),]
rownames(ref_G)=ref_G[,1]

Specify Contrast(s) of interest

Use a linear model to fit voomed counts (logCPM). Instanciate the comparisons of interest by using the contrast approach from limma. Extract the results, annotate the tables and merge them into a global result dataframe.

Comparisons between conditions Batch3

fitB3 <- lmFit(yB3, mmB3)
## Coefficients not estimable: sampleAB253 sampleAB257
contr <- makeContrasts(
  A = lesion_typeWML - lesion_typeNAWM,
  B = lesion_typeWML - lesion_typeCL,
  C = lesion_typeCL - lesion_typeNAGM,
  D = lesion_typeNAWM - lesion_typeNAGM,
  E = lesion_typePVML - lesion_typeNAWM,
  F = lesion_typePVML - lesion_typeWML,
  G = lesion_typePVML - lesion_typeCL,
  H = lesion_typeNAGM - lesion_typeGM,
  I = lesion_typePVML - lesion_typePVM,
  J = lesion_typeWML - lesion_typeWM,
  levels = colnames(coef(fitB3)))
tmp <- contrasts.fit(fitB3, contr)
tmp <- eBayes(tmp)

resA_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='A')
resB_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='B')
resC_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='C')
resD_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='D')
resE_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='E')
resF_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='F')
resG_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='G')
resH_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='H')
resI_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='I')
resJ_B3 <- topTable(tmp, sort.by = "P", n = Inf,coef='J')

Comparisons between conditions All data

fit <- lmFit(y, mm)
## Coefficients not estimable: sampleAB257 sexM
contr <- makeContrasts(
  A = lesion_typeWML - lesion_typeNAWM,
  B = lesion_typeWML - lesion_typeCL,
  C = lesion_typeCL - lesion_typeNAGM,
  D = lesion_typeNAWM - lesion_typeNAGM,
  E = lesion_typePVML - lesion_typeNAWM,
  F = lesion_typePVML - lesion_typeWML,
  G = lesion_typePVML - lesion_typeCL,
  H = lesion_typeNAGM - lesion_typeGM,
  I = lesion_typePVML - lesion_typePVM,
  J = lesion_typeWML - lesion_typeWM,
  levels = colnames(coef(fit)))
tmp <- contrasts.fit(fit, contr)
tmp <- eBayes(tmp)

resA <- topTable(tmp, sort.by = "P", n = Inf,coef='A')
resB <- topTable(tmp, sort.by = "P", n = Inf,coef='B')
resC <- topTable(tmp, sort.by = "P", n = Inf,coef='C')
resD <- topTable(tmp, sort.by = "P", n = Inf,coef='D')
resE <- topTable(tmp, sort.by = "P", n = Inf,coef='E')
resF <- topTable(tmp, sort.by = "P", n = Inf,coef='F')
resG <- topTable(tmp, sort.by = "P", n = Inf,coef='G')
resH <- topTable(tmp, sort.by = "P", n = Inf,coef='H')
resI <- topTable(tmp, sort.by = "P", n = Inf,coef='I')
resJ <- topTable(tmp, sort.by = "P", n = Inf,coef='J')

Comparisons between group MS & HC

mmG <- model.matrix(~ 0 + group+batch+sample+sex, data=meta[colnames(d),])

# Filtering low expressed genes
yG <- voom(d[,rownames(meta)], mm, plot = F)
## Coefficients not estimable: sampleAB257 sexM
fitG <- lmFit(yG, mmG)
## Coefficients not estimable: sampleAB257 sexM
contrG <- makeContrasts(
  A = groupMS.WM - groupHC.WM,
  B = groupMS.Cortical - groupHC.Cortical,
  C = groupMS.PV - groupHC.PV,
  levels = colnames(coef(fitG)))
tmpG <- contrasts.fit(fitG, contrG)
tmpG <- eBayes(tmpG)

resA_group <- topTable(tmpG, sort.by = "P", n = Inf,coef='A')
resB_group <- topTable(tmpG, sort.by = "P", n = Inf,coef='B')
resC_group <- topTable(tmpG, sort.by = "P", n = Inf,coef='C')

Batch3

Comparisons between conditions Batch3 (Volcano plot)

Plot results in the form of volcano plots. Thresholds used are adj.P.Val<0.05 & |logFC|>1. Red genes are upregulated in the first group as found in the subtitle of each graph.

PVML_vs_NAWM

resE_B3$gene_id =  rownames(resE_B3)
resE_B3$comp = 'PVML_vs_NAWM'
resE_B3$gene_name = ref_G[rownames(resE_B3),3]
resE_B3=resE_B3[order(-abs(resE_B3$logFC)),]
resE_B3$show = F
resE_B3[1:15,'show']=T
resE_B3$sign = 'NS'
resE_B3[which(resE_B3$adj.P.Val<0.05&resE_B3$logFC>1),'sign']='UP'
resE_B3[which(resE_B3$adj.P.Val<0.05&resE_B3$logFC<(-1)),'sign']='DN'
ggplot(resE_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resE_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_NAWM')

WML_vs_NAWM

resA_B3$gene_id =  rownames(resA_B3)
resA_B3$comp = 'WML_vs_NAWM'
resA_B3$gene_name = ref_G[rownames(resA_B3),3]
resA_B3=resA_B3[order(-abs(resA_B3$logFC)),]
resA_B3$show = F
resA_B3[1:15,'show']=T
resA_B3$sign = 'NS'
resA_B3[which(resA_B3$adj.P.Val<0.05&resA_B3$logFC>1),'sign']='UP'
resA_B3[which(resA_B3$adj.P.Val<0.05&resA_B3$logFC<(-1)),'sign']='DN'
ggplot(resA_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resA_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('WML_vs_NAWM')

WML_vs_CL

resB_B3$gene_id =  rownames(resB_B3)
resB_B3$comp = 'WML_vs_CL'
resB_B3$gene_name = ref_G[rownames(resB_B3),3]
resB_B3=resB_B3[order(-abs(resB_B3$logFC)),]
resB_B3$show = F
resB_B3[1:15,'show']=T
resB_B3$sign = 'NS'
resB_B3[which(resB_B3$adj.P.Val<0.05&resB_B3$logFC>1),'sign']='UP'
resB_B3[which(resB_B3$adj.P.Val<0.05&resB_B3$logFC<(-1)),'sign']='DN'
ggplot(resB_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resB_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('WML_vs_CL')

CL_vs_NAGM

resC_B3$gene_id =  rownames(resC_B3)
resC_B3$comp = 'CL_vs_NAGM'
resC_B3$gene_name = ref_G[rownames(resC_B3),3]
resC_B3=resC_B3[order(-abs(resC_B3$logFC)),]
resC_B3$show = F
resC_B3[1:15,'show']=T
resC_B3$sign = 'NS'
resC_B3[which(resC_B3$adj.P.Val<0.05&resC_B3$logFC>1),'sign']='UP'
resC_B3[which(resC_B3$adj.P.Val<0.05&resC_B3$logFC<(-1)),'sign']='DN'
ggplot(resC_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resC_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('CL_vs_NAGM')

NAWM_vs_NAGM

resD_B3$gene_id =  rownames(resD_B3)
resD_B3$comp = 'NAWM_vs_NAGM'
resD_B3$gene_name = ref_G[rownames(resD_B3),3]
resD_B3=resD_B3[order(-abs(resD_B3$logFC)),]
resD_B3$show = F
resD_B3[1:15,'show']=T
resD_B3$sign = 'NS'
resD_B3[which(resD_B3$adj.P.Val<0.05&resD_B3$logFC>1),'sign']='UP'
resD_B3[which(resD_B3$adj.P.Val<0.05&resD_B3$logFC<(-1)),'sign']='DN'
ggplot(resD_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resD_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('NAWM_vs_NAGM')

PVML_vs_WML

resF_B3$gene_id =  rownames(resF_B3)
resF_B3$comp = 'PVML_vs_WML'
resF_B3$gene_name = ref_G[rownames(resF_B3),3]
resF_B3=resF_B3[order(-abs(resF_B3$logFC)),]
resF_B3$show = F
resF_B3[1:15,'show']=T
resF_B3$sign = 'NS'
resF_B3[which(resF_B3$adj.P.Val<0.05&resF_B3$logFC>1),'sign']='UP'
resF_B3[which(resF_B3$adj.P.Val<0.05&resF_B3$logFC<(-1)),'sign']='DN'
ggplot(resF_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resF_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_WML')

PVML_vs_CL

resG_B3$gene_id =  rownames(resG_B3)
resG_B3$comp = 'PVML_vs_CL'
resG_B3$gene_name = ref_G[rownames(resG_B3),3]
resG_B3=resG_B3[order(-abs(resG_B3$logFC)),]
resG_B3$show = F
resG_B3[1:15,'show']=T
resG_B3$sign = 'NS'
resG_B3[which(resG_B3$adj.P.Val<0.05&resG_B3$logFC>1),'sign']='UP'
resG_B3[which(resG_B3$adj.P.Val<0.05&resG_B3$logFC<(-1)),'sign']='DN'
ggplot(resG_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resG_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_CL')

NAGM_vs_GM

resH_B3$gene_id =  rownames(resH_B3)
resH_B3$comp = 'NAGM_vs_GM'
resH_B3$gene_name = ref_G[rownames(resH_B3),3]
resH_B3=resH_B3[order(-abs(resH_B3$logFC)),]
resH_B3$show = F
resH_B3[1:15,'show']=T
resH_B3$sign = 'NS'
resH_B3[which(resH_B3$adj.P.Val<0.05&resH_B3$logFC>1),'sign']='UP'
resH_B3[which(resH_B3$adj.P.Val<0.05&resH_B3$logFC<(-1)),'sign']='DN'
ggplot(resH_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resH_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('NAGM_vs_GM')

PVML_vs_PVM

resI_B3$gene_id =  rownames(resI_B3)
resI_B3$comp = 'PVML_vs_PVM'
resI_B3$gene_name = ref_G[rownames(resI_B3),3]
resI_B3=resI_B3[order(-abs(resI_B3$logFC)),]
resI_B3$show = F
resI_B3[1:15,'show']=T
resI_B3$sign = 'NS'
resI_B3[which(resI_B3$adj.P.Val<0.05&resI_B3$logFC>1),'sign']='UP'
resI_B3[which(resI_B3$adj.P.Val<0.05&resI_B3$logFC<(-1)),'sign']='DN'
ggplot(resI_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resI_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_PVM')

WML_vs_WM

resJ_B3$gene_id =  rownames(resJ_B3)
resJ_B3$comp = 'WML_vs_WM'
resJ_B3$gene_name = ref_G[rownames(resJ_B3),3]
resJ_B3=resJ_B3[order(-abs(resJ_B3$logFC)),]
resJ_B3$show = F
resJ_B3[1:15,'show']=T
resJ_B3$sign = 'NS'
resJ_B3[which(resJ_B3$adj.P.Val<0.05&resJ_B3$logFC>1),'sign']='UP'
resJ_B3[which(resJ_B3$adj.P.Val<0.05&resJ_B3$logFC<(-1)),'sign']='DN'
ggplot(resJ_B3,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resJ_B3,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('WML_vs_WM')

ALL DATA

Comparisons between conditions (Volcano plot)

Plot results in the form of volcano plots. Thresholds used are adj.P.Val<0.05 & |logFC|>1. Red genes are upregulated in the first group as found in the subtitle of each graph.

PVML_vs_NAWM

resE$gene_id =  rownames(resE)
resE$comp = 'PVML_vs_NAWM'
resE$gene_name = ref_G[rownames(resE),3]
resE=resE[order(-abs(resE$logFC)),]
resE$show = F
resE[1:15,'show']=T
resE$sign = 'NS'
resE[which(resE$adj.P.Val<0.05&resE$logFC>1),'sign']='UP'
resE[which(resE$adj.P.Val<0.05&resE$logFC<(-1)),'sign']='DN'
ggplot(resE,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resE,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_NAWM')

WML_vs_NAWM

resA$gene_id =  rownames(resA)
resA$comp = 'WML_vs_NAWM'
resA$gene_name = ref_G[rownames(resA),3]
resA=resA[order(-abs(resA$logFC)),]
resA$show = F
resA[1:15,'show']=T
resA$sign = 'NS'
resA[which(resA$adj.P.Val<0.05&resA$logFC>1),'sign']='UP'
resA[which(resA$adj.P.Val<0.05&resA$logFC<(-1)),'sign']='DN'
ggplot(resA,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resA,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('WML_vs_NAWM')

WML_vs_CL

resB$gene_id =  rownames(resB)
resB$comp = 'WML_vs_CL'
resB$gene_name = ref_G[rownames(resB),3]
resB=resB[order(-abs(resB$logFC)),]
resB$show = F
resB[1:15,'show']=T
resB$sign = 'NS'
resB[which(resB$adj.P.Val<0.05&resB$logFC>1),'sign']='UP'
resB[which(resB$adj.P.Val<0.05&resB$logFC<(-1)),'sign']='DN'
ggplot(resB,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resB,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('WML_vs_CL')

CL_vs_NAGM

resC$gene_id =  rownames(resC)
resC$comp = 'CL_vs_NAGM'
resC$gene_name = ref_G[rownames(resC),3]
resC=resC[order(-abs(resC$logFC)),]
resC$show = F
resC[1:15,'show']=T
resC$sign = 'NS'
resC[which(resC$adj.P.Val<0.05&resC$logFC>1),'sign']='UP'
resC[which(resC$adj.P.Val<0.05&resC$logFC<(-1)),'sign']='DN'
ggplot(resC,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resC,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('CL_vs_NAGM')

NAWM_vs_NAGM

resD$gene_id =  rownames(resD)
resD$comp = 'NAWM_vs_NAGM'
resD$gene_name = ref_G[rownames(resD),3]
resD=resD[order(-abs(resD$logFC)),]
resD$show = F
resD[1:15,'show']=T
resD$sign = 'NS'
resD[which(resD$adj.P.Val<0.05&resD$logFC>1),'sign']='UP'
resD[which(resD$adj.P.Val<0.05&resD$logFC<(-1)),'sign']='DN'
ggplot(resD,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resD,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('NAWM_vs_NAGM')

PVML_vs_WML

resF$gene_id =  rownames(resF)
resF$comp = 'PVML_vs_WML'
resF$gene_name = ref_G[rownames(resF),3]
resF=resF[order(-abs(resF$logFC)),]
resF$show = F
resF[1:15,'show']=T
resF$sign = 'NS'
resF[which(resF$adj.P.Val<0.05&resF$logFC>1),'sign']='UP'
resF[which(resF$adj.P.Val<0.05&resF$logFC<(-1)),'sign']='DN'
ggplot(resF,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resF,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_WML')

PVML_vs_CL

resG$gene_id =  rownames(resG)
resG$comp = 'PVML_vs_CL'
resG$gene_name = ref_G[rownames(resG),3]
resG=resG[order(-abs(resG$logFC)),]
resG$show = F
resG[1:15,'show']=T
resG$sign = 'NS'
resG[which(resG$adj.P.Val<0.05&resG$logFC>1),'sign']='UP'
resG[which(resG$adj.P.Val<0.05&resG$logFC<(-1)),'sign']='DN'
ggplot(resG,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resG,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_CL')

NAGM_vs_GM

resH$gene_id =  rownames(resH)
resH$comp = 'NAGM_vs_GM'
resH$gene_name = ref_G[rownames(resH),3]
resH=resH[order(-abs(resH$logFC)),]
resH$show = F
resH[1:15,'show']=T
resH$sign = 'NS'
resH[which(resH$adj.P.Val<0.05&resH$logFC>1),'sign']='UP'
resH[which(resH$adj.P.Val<0.05&resH$logFC<(-1)),'sign']='DN'
ggplot(resH,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resH,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('NAGM_vs_GM')

PVML_vs_PVM

resI$gene_id =  rownames(resI)
resI$comp = 'PVML_vs_PVM'
resI$gene_name = ref_G[rownames(resI),3]
resI=resI[order(-abs(resI$logFC)),]
resI$show = F
resI[1:15,'show']=T
resI$sign = 'NS'
resI[which(resI$adj.P.Val<0.05&resI$logFC>1),'sign']='UP'
resI[which(resI$adj.P.Val<0.05&resI$logFC<(-1)),'sign']='DN'
ggplot(resI,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resI,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('PVML_vs_PVM')

WML_vs_WM

resJ$gene_id =  rownames(resJ)
resJ$comp = 'WML_vs_WM'
resJ$gene_name = ref_G[rownames(resJ),3]
resJ=resJ[order(-abs(resJ$logFC)),]
resJ$show = F
resJ[1:15,'show']=T
resJ$sign = 'NS'
resJ[which(resJ$adj.P.Val<0.05&resJ$logFC>1),'sign']='UP'
resJ[which(resJ$adj.P.Val<0.05&resJ$logFC<(-1)),'sign']='DN'
ggplot(resJ,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resJ,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('WML_vs_WM')

Comparisons between group MS & HC (Volcano plot)

Plot results in the form of volcano plots. Thresholds used are adj.P.Val<0.05 & |logFC|>1. Red genes are upregulated in the first group as found in the subtitle of each graph.

groupMS.WM_vs_groupHC.WM

resA_group$gene_id =  rownames(resA_group)
resA_group$comp = 'groupMS.WM_vs_groupHC.WM'
resA_group$gene_name = ref_G[rownames(resA_group),3]
resA_group=resA_group[order(-abs(resA_group$logFC)),]
resA_group$show = F
resA_group[1:15,'show']=T
resA_group$sign = 'NS'
resA_group[which(resA_group$adj.P.Val<0.05&resA_group$logFC>1),'sign']='UP'
resA_group[which(resA_group$adj.P.Val<0.05&resA_group$logFC<(-1)),'sign']='DN'
ggplot(resA_group,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resA_group,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('groupMS.WM_vs_groupHC.WM')

groupMS.Cortical_vs_groupHC.Cortical

resB_group$gene_id =  rownames(resB_group)
resB_group$comp = 'groupMS.Cortical_vs_groupHC.Cortical'
resB_group$gene_name = ref_G[rownames(resB_group),3]
resB_group=resB_group[order(-abs(resB_group$logFC)),]
resB_group$show = F
resB_group[1:15,'show']=T
resB_group$sign = 'NS'
resB_group[which(resB_group$adj.P.Val<0.05&resB_group$logFC>1),'sign']='UP'
resB_group[which(resB_group$adj.P.Val<0.05&resB_group$logFC<(-1)),'sign']='DN'
ggplot(resB_group,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resB_group,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('groupMS.Cortical_vs_groupHC.Cortical')

groupMS.PV_vs_groupHC.PV

resC_group$gene_id =  rownames(resC_group)
resC_group$comp = 'groupMS.PV_vs_groupHC.PV'
resC_group$gene_name = ref_G[rownames(resC_group),3]
resC_group=resC_group[order(-abs(resC_group$logFC)),]
resC_group$show = F
resC_group[1:15,'show']=T
resC_group$sign = 'NS'
resC_group[which(resC_group$adj.P.Val<0.05&resC_group$logFC>1),'sign']='UP'
resC_group[which(resC_group$adj.P.Val<0.05&resC_group$logFC<(-1)),'sign']='DN'
ggplot(resC_group,aes(logFC,-log10(P.Value),col=sign))+geom_point_interactive(aes(x=logFC,y=-log10(P.Value),color=sign,tooltip = gene_name,data_id = gene_name))+theme_bw()+scale_color_manual(values=sign_cols)+ geom_label_repel(data=subset(resC_group,show==T),aes(label=gene_name),col='black',size=3,force=30,fill = alpha(c("white"),0.5))+ggtitle('groupMS.PV_vs_groupHC.PV')

resE100=head(resE,100)
pheatmap::pheatmap(y$E[rownames(resE100),],scale = "row")

pheatmap::pheatmap(cor(y$E[rownames(resE100),]))

Gene ontology term enrichment analysis of comparisons between conditions

This extracts biological signatures from the transcriptomic changes found in the differential expression analysis

enr.a=find_enrichments(resA)
enr.b=find_enrichments(resB)
enr.c=find_enrichments(resC)
enr.d=find_enrichments(resD)
enr.e=find_enrichments(resE)
enr.f=find_enrichments(resF)
enr.g=find_enrichments(resG)
enr.h=find_enrichments(resH)
enr.i=find_enrichments(resI)
enr.j=find_enrichments(resJ)

Biological Processes

enrA = enr.a$enr$bp
enrA=enrA[order(enrA$padj),]
enrA$comp = unique(resA$comp)

enrB = enr.b$enr$bp
enrB=enrB[order(enrB$padj),]
enrB$comp = unique(resB$comp)

enrC = enr.c$enr$bp
enrC=enrC[order(enrC$padj),]
enrC$comp = unique(resC$comp)

enrD = enr.d$enr$bp
enrD=enrD[order(enrD$padj),]
enrD$comp = unique(resD$comp)

enrE = enr.e$enr$bp
enrE=enrE[order(enrE$padj),]
enrE$comp = unique(resE$comp)

enrF = enr.f$enr$bp
enrF=enrF[order(enrF$padj),]
enrF$comp = unique(resF$comp)

enrG = enr.g$enr$bp
enrG=enrG[order(enrG$padj),]
enrG$comp = unique(resG$comp)

enrH = enr.h$enr$bp
enrH=enrH[order(enrH$padj),]
enrH$comp = unique(resH$comp)

enrI = enr.i$enr$bp
enrI=enrI[order(enrI$padj),]
enrI$comp = unique(resI$comp)

enrJ = enr.j$enr$bp
enrJ=enrJ[order(enrJ$padj),]
enrJ$comp = unique(resJ$comp)

paths = c(head(enrA,10)$pathway,
          head(enrB,10)$pathway,
          head(enrC,10)$pathway,
          head(enrD,10)$pathway,
          head(enrE,10)$pathway,
          head(enrF,10)$pathway,
          head(enrG,10)$pathway,
          head(enrH,10)$pathway,
          head(enrI,10)$pathway,
          head(enrJ,10)$pathway)

paths=unique(paths)
rb=rbind(enrA[which(enrA$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrB[which(enrB$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrC[which(enrC$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrD[which(enrD$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrE[which(enrE$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrF[ which(enrF$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrG[which(enrG$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrH[which(enrH$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrI[ which(enrI$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrJ[which(enrJ$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')])
rb[which(rb$padj>=0.1),'NES']=NA
rb$pathway <- gsub("GOBP", "",rb$pathway)
rb$pathway <- gsub("_", " ",rb$pathway)
ss  = strsplit(as.character(rb$leadingEdge),',')
rb$n_genes = unlist(lapply(ss,length))
rb$pval <- as.numeric(rb$pval)
rb[which(rb$pval>=0.1),'pval']=NA
ggplot(na.omit(rb),aes(pathway,comp,fill=NES))+geom_point(aes(size=pval),shape=21)+scale_x_discrete(label = function(x) stringr::str_trunc(x,50))+theme_bw()+coord_flip()+theme(axis.text.x=element_text( angle=60,hjust=1))+scale_fill_gradientn(colors=rev(brewer.pal(11,'RdYlBu')))+ scale_size(trans = 'reverse')+ theme(axis.text = element_text(face="bold"))

Molecular Function

enrA = enr.a$enr$mf
enrA=enrA[order(enrA$padj),]
enrA$comp = unique(resA$comp)

enrB = enr.b$enr$mf
enrB=enrB[order(enrB$padj),]
enrB$comp = unique(resB$comp)

enrC = enr.c$enr$mf
enrC=enrC[order(enrC$padj),]
enrC$comp = unique(resC$comp)

enrD = enr.d$enr$mf
enrD=enrD[order(enrD$padj),]
enrD$comp = unique(resD$comp)

enrE = enr.e$enr$mf
enrE=enrE[order(enrE$padj),]
enrE$comp = unique(resE$comp)

enrF = enr.f$enr$mf
enrF=enrF[order(enrF$padj),]
enrF$comp = unique(resF$comp)

enrG = enr.g$enr$mf
enrG=enrG[order(enrG$padj),]
enrG$comp = unique(resG$comp)

enrH = enr.h$enr$mf
enrH=enrH[order(enrH$padj),]
enrH$comp = unique(resH$comp)

enrI = enr.i$enr$mf
enrI=enrI[order(enrI$padj),]
enrI$comp = unique(resI$comp)

enrJ = enr.j$enr$mf
enrJ=enrJ[order(enrJ$padj),]
enrJ$comp = unique(resJ$comp)

paths = c(head(enrA,10)$pathway,
          head(enrB,10)$pathway,
          head(enrC,10)$pathway,
          head(enrD,10)$pathway,
          head(enrE,10)$pathway,
          head(enrF,10)$pathway,
          head(enrG,10)$pathway,
          head(enrH,10)$pathway,
          head(enrI,10)$pathway,
          head(enrJ,10)$pathway)

paths=unique(paths)
rb=rbind(enrA[which(enrA$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrB[which(enrB$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrC[which(enrC$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrD[which(enrD$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrE[which(enrE$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrF[ which(enrF$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrG[which(enrG$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrH[which(enrH$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrI[ which(enrI$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrJ[which(enrJ$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')])

rb[which(rb$padj>=0.1),'NES']=NA
rb$pathway <- gsub("GOMF", "",rb$pathway)
rb$pathway <- gsub("_", " ",rb$pathway)
ss  = strsplit(as.character(rb$leadingEdge),',')
rb$n_genes = unlist(lapply(ss,length))
rb$pval <- as.numeric(rb$pval)
rb[which(rb$pval>=0.1),'pval']=NA
ggplot(na.omit(rb),aes(pathway,comp,fill=NES))+geom_point(aes(size=pval),shape=21)+scale_x_discrete(label = function(x) stringr::str_trunc(x,50))+theme_bw()+coord_flip()+theme(axis.text.x=element_text( angle=60,hjust=1))+scale_fill_gradientn(colors=rev(brewer.pal(11,'RdYlBu')))+ scale_size(trans = 'reverse')+ theme(axis.text = element_text(face="bold"))

KEGG pathway database

enrA = enr.a$enr$kegg
enrA=enrA[order(enrA$padj),]
enrA$comp = unique(resA$comp)

enrB = enr.b$enr$kegg
enrB=enrB[order(enrB$padj),]
enrB$comp = unique(resB$comp)

enrC = enr.c$enr$kegg
enrC=enrC[order(enrC$padj),]
enrC$comp = unique(resC$comp)

enrD = enr.d$enr$kegg
enrD=enrD[order(enrD$padj),]
enrD$comp = unique(resD$comp)

enrE = enr.e$enr$kegg
enrE=enrE[order(enrE$padj),]
enrE$comp = unique(resE$comp)

enrF = enr.f$enr$kegg
enrF=enrF[order(enrF$padj),]
enrF$comp = unique(resF$comp)

enrG = enr.g$enr$kegg
enrG=enrG[order(enrG$padj),]
enrG$comp = unique(resG$comp)

enrH = enr.h$enr$kegg
enrH=enrH[order(enrH$padj),]
enrH$comp = unique(resH$comp)

enrI = enr.i$enr$kegg
enrI=enrI[order(enrI$padj),]
enrI$comp = unique(resI$comp)

enrJ = enr.j$enr$kegg
enrJ=enrJ[order(enrJ$padj),]
enrJ$comp = unique(resJ$comp)

paths = c(head(enrA,10)$pathway,
          head(enrB,10)$pathway,
          head(enrC,10)$pathway,
          head(enrD,10)$pathway,
          head(enrE,10)$pathway,
          head(enrF,10)$pathway,
          head(enrG,10)$pathway,
          head(enrH,10)$pathway,
          head(enrI,10)$pathway,
          head(enrJ,10)$pathway)

paths=unique(paths)
rb=rbind(enrA[which(enrA$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrB[which(enrB$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrC[which(enrC$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrD[which(enrD$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrE[which(enrE$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrF[ which(enrF$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrG[which(enrG$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrH[which(enrH$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrI[ which(enrI$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')],
         enrJ[which(enrJ$pathway%in%paths),c('pval','leadingEdge','NES','comp','pathway')])

rb[which(rb$padj>=0.1),'NES']=NA
rb$pathway <- gsub("KEGG", "",rb$pathway)
rb$pathway <- gsub("_", " ",rb$pathway)
ss  = strsplit(as.character(rb$leadingEdge),',')
rb$n_genes = unlist(lapply(ss,length))
rb$pval <- as.numeric(rb$pval)
rb[which(rb$pval>=0.1),'pval']=NA
ggplot(na.omit(rb),aes(pathway,comp,fill=NES))+geom_point(aes(size=pval),shape=21)+scale_x_discrete(label = function(x) stringr::str_trunc(x,50))+theme_bw()+coord_flip()+theme(axis.text.x=element_text( angle=60,hjust=1))+scale_fill_gradientn(colors=rev(brewer.pal(11,'RdYlBu')))+ scale_size(trans = 'reverse')+ theme(axis.text = element_text(face="bold"))

Plot the overlap between DEGs with an upset plot/Venn Diagram

Overlap of DEGs in comparaisons

res=rbind(resA,resB,resC,resD,resE,resF,resG,resH,resI,resJ)
res= na.omit(res)

degsListvs = list(
    WML_vs_NAWM=subset(res,comp=='WML_vs_NAWM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    PVML_vs_WML=subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    PVML_vs_NAWM=subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    WML_vs_CL=subset(res,comp=='WML_vs_CL' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    #NAWM_vs_NAGM=subset(res,comp=='NAWM_vs_NAGM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    NAGM_vs_GM=subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    PVML_vs_PVM=subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    WML_vs_WM=subset(res,comp=='WML_vs_WM' & P.Value<0.05 & abs(logFC)>1)$gene_name)
    

degsListUP = list(
    WML_vs_NAWM=subset(res,comp=='WML_vs_NAWM' & P.Value<0.05 & (logFC)>1)$gene_name,
    PVML_vs_WML=subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & (logFC)>1)$gene_name,
    PVML_vs_NAWM=subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & (logFC)>1)$gene_name,
    WML_vs_CL=subset(res,comp=='WML_vs_CL' & P.Value<0.05 & (logFC)>1)$gene_name,
    #NAWM_vs_NAGM=subset(res,comp=='NAWM_vs_NAGM' & P.Value<0.05 & (logFC)>1)$gene_name,
    NAGM_vs_GM=subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & (logFC)>1)$gene_name,
    PVML_vs_PVM=subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & (logFC)>1)$gene_name,
    WML_vs_WM=subset(res,comp=='WML_vs_WM' & P.Value<0.05 & (logFC)>1)$gene_name)

degsListDOWN = list(
    WML_vs_NAWM=subset(res,comp=='WML_vs_NAWM' & P.Value<0.05 & (logFC)<1)$gene_name,
    PVML_vs_WML=subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & (logFC)<1)$gene_name,
    PVML_vs_NAWM=subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & (logFC)<1)$gene_name,
    WML_vs_CL=subset(res,comp=='WML_vs_CL' & P.Value<0.05 & (logFC)<1)$gene_name,
    #NAWM_vs_NAGM=subset(res,comp=='NAWM_vs_NAGM' & P.Value<0.05 & (logFC)<1)$gene_name,
    NAGM_vs_GM=subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & (logFC)<1)$gene_name,
    PVML_vs_PVM=subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & (logFC)<1)$gene_name,
    WML_vs_WM=subset(res,comp=='WML_vs_WM' & P.Value<0.05 & (logFC)<1)$gene_name)

degsListvs_UpSet = list(
    WML_vs_NAWM=subset(res,comp=='WML_vs_NAWM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    PVML_vs_WML=subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    PVML_vs_NAWM=subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    WML_vs_CL=subset(res,comp=='WML_vs_CL' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    NAWM_vs_NAGM=subset(res,comp=='NAWM_vs_NAGM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    NAGM_vs_GM=subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    PVML_vs_PVM=subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & abs(logFC)>1)$gene_name,
    WML_vs_WM=subset(res,comp=='WML_vs_WM' & P.Value<0.05 & abs(logFC)>1)$gene_name)
    

degsListUP_UpSet = list(
    WML_vs_NAWM=subset(res,comp=='WML_vs_NAWM' & P.Value<0.05 & (logFC)>1)$gene_name,
    PVML_vs_WML=subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & (logFC)>1)$gene_name,
    PVML_vs_NAWM=subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & (logFC)>1)$gene_name,
    WML_vs_CL=subset(res,comp=='WML_vs_CL' & P.Value<0.05 & (logFC)>1)$gene_name,
    NAWM_vs_NAGM=subset(res,comp=='NAWM_vs_NAGM' & P.Value<0.05 & (logFC)>1)$gene_name,
    NAGM_vs_GM=subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & (logFC)>1)$gene_name,
    PVML_vs_PVM=subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & (logFC)>1)$gene_name,
    WML_vs_WM=subset(res,comp=='WML_vs_WM' & P.Value<0.05 & (logFC)>1)$gene_name)

degsListDOWN_UpSet = list(
    WML_vs_NAWM=subset(res,comp=='WML_vs_NAWM' & P.Value<0.05 & (logFC)<1)$gene_name,
    PVML_vs_WML=subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & (logFC)<1)$gene_name,
    PVML_vs_NAWM=subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & (logFC)<1)$gene_name,
    WML_vs_CL=subset(res,comp=='WML_vs_CL' & P.Value<0.05 & (logFC)<1)$gene_name,
    NAWM_vs_NAGM=subset(res,comp=='NAWM_vs_NAGM' & P.Value<0.05 & (logFC)<1)$gene_name,
    NAGM_vs_GM=subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & (logFC)<1)$gene_name,
    PVML_vs_PVM=subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & (logFC)<1)$gene_name,
    WML_vs_WM=subset(res,comp=='WML_vs_WM' & P.Value<0.05 & (logFC)<1)$gene_name)

Venn Diagram

Overlap between DEGs (All)

#label = "count"
ggVennDiagram(degsListvs,label = "percent", label_alpha = 0)+ scale_fill_distiller(palette = "Reds", direction = 1)

Overlap between DEGs (UP)

#label = "count"
ggVennDiagram(degsListUP,label = "percent", label_alpha = 0)+ scale_fill_distiller(palette = "Reds", direction = 1)

Overlap between DEGs (DOWN)

#label = "count"
ggVennDiagram(degsListDOWN,label = "percent", label_alpha = 0)+ scale_fill_distiller(palette = "Reds", direction = 1)

UpSet plot

Overlap between DEGs (All)

upset(fromList(degsListvs_UpSet),sets = c('WML_vs_NAWM','PVML_vs_WML','PVML_vs_NAWM',"WML_vs_CL","NAWM_vs_NAGM","NAGM_vs_GM", "PVML_vs_PVM", "WML_vs_WM"),order.by='freq',decreasing=c(F), text.scale = c(1.3, 2, 1.5, 1.5, 2, 2))

Overlap between DEGs (UP)

upset(fromList(degsListUP_UpSet),sets = c('WML_vs_NAWM','PVML_vs_WML','PVML_vs_NAWM',"WML_vs_CL","NAWM_vs_NAGM","NAGM_vs_GM", "PVML_vs_PVM", "WML_vs_WM"),order.by='freq',decreasing = c(F),text.scale = c(1.3, 2, 1.5, 1.5, 2, 2))

Overlap between DEGs (DOWN)

upset(fromList(degsListDOWN_UpSet),sets = c('WML_vs_NAWM','PVML_vs_WML','PVML_vs_NAWM',"WML_vs_CL","NAWM_vs_NAGM","NAGM_vs_GM", "PVML_vs_PVM", "WML_vs_WM"),order.by='freq',decreasing = c(F),text.scale = c(1.3, 2, 1.5, 1.5, 2, 2))

Gene Ontology enrichment analysis of overlapping genes in comparisons between conditions(The top 5 specific Genesets from the upsets plot)

WML_vs_CL=na.omit(subset(res,comp=='WML_vs_CL' & P.Value<0.05 & logFC>1)$gene_name)
PVML_vs_PVM=na.omit(subset(res,comp=='PVML_vs_PVM' & P.Value<0.05 & logFC>1)$gene_name)
PVML_vs_NAWM=na.omit(subset(res,comp=='PVML_vs_NAWM' & P.Value<0.05 & logFC>1)$gene_name)
NAGM_vs_GM=na.omit(subset(res,comp=='NAGM_vs_GM' & P.Value<0.05 & logFC>1)$gene_name)
PVML_vs_WML=na.omit(subset(res,comp=='PVML_vs_WML' & P.Value<0.05 & logFC>1)$gene_name)

specificWML_vs_CL = as.character(WML_vs_CL[!WML_vs_CL%in%PVML_vs_PVM & !WML_vs_CL%in%PVML_vs_NAWM & !WML_vs_CL%in%NAGM_vs_GM & !WML_vs_CL%in%PVML_vs_WML])

specificPVML_vs_PVM = as.character(PVML_vs_PVM[!PVML_vs_PVM%in%WML_vs_CL & !PVML_vs_PVM%in%PVML_vs_NAWM & !PVML_vs_PVM%in%NAGM_vs_GM & !PVML_vs_PVM%in%PVML_vs_WML])

specificPVML_vs_NAWM = as.character(PVML_vs_NAWM[!PVML_vs_NAWM%in%WML_vs_CL & !PVML_vs_NAWM%in%PVML_vs_PVM & !PVML_vs_NAWM%in%NAGM_vs_GM & !PVML_vs_NAWM%in%PVML_vs_WML])

specificNAGM_vs_GM = as.character(NAGM_vs_GM[!NAGM_vs_GM%in%WML_vs_CL & !NAGM_vs_GM%in%PVML_vs_PVM & !NAGM_vs_GM%in%PVML_vs_NAWM & !NAGM_vs_GM%in%PVML_vs_WML])
  
commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = as.character(PVML_vs_NAWM[!PVML_vs_NAWM%in%WML_vs_CL & PVML_vs_NAWM%in%PVML_vs_WML & PVML_vs_NAWM%in%PVML_vs_PVM & !PVML_vs_NAWM%in%NAGM_vs_GM])
  

list = list(commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM=commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM, specificWML_vs_CL=specificWML_vs_CL,specificPVML_vs_PVM=specificPVML_vs_PVM, specificPVML_vs_NAWM=specificPVML_vs_NAWM,specificNAGM_vs_GM=specificNAGM_vs_GM)

fishBP = fisher_enrichment(list,go,unique(ref_G$SYMBOL))
fishMF = fisher_enrichment(list,mf,unique(ref_G$SYMBOL))
#fishKEGG = fisher_enrichment(list,kegg,unique(ref_G$SYMBOL))

biological processes

  enr.specificWML_vs_CL = fishBP$specificWML_vs_CL
  enr.specificWML_vs_CL$pathway = rownames(enr.specificWML_vs_CL)
  enr.specificWML_vs_CL$fdr <- as.numeric(enr.specificWML_vs_CL$fdr) 
  enr.specificWML_vs_CL$OR.odds.ratio <- as.numeric(enr.specificWML_vs_CL$OR.odds.ratio) 
  enr.specificWML_vs_CL$NOM_pval <- as.numeric(enr.specificWML_vs_CL$NOM_pval)
  topspecificWML_vs_CL <- topPath(fishBP ,"specificWML_vs_CL", 0.5, 20)
  enr.specificWML_vs_CL <- enr.specificWML_vs_CL[match(topspecificWML_vs_CL, rownames(enr.specificWML_vs_CL)),]
  rownames(enr.specificWML_vs_CL) <- gsub("GOBP_", "",rownames(enr.specificWML_vs_CL)) 
  enr.specificWML_vs_CL=enr.specificWML_vs_CL[order(enr.specificWML_vs_CL$NOM_pval),]
  enr.specificWML_vs_CL$comp = "specificWML_vs_CL"
  
  enr.specificPVML_vs_PVM = fishBP$specificPVML_vs_PVM
  enr.specificPVML_vs_PVM$pathway = rownames(enr.specificPVML_vs_PVM)
  enr.specificPVML_vs_PVM$fdr <- as.numeric(enr.specificPVML_vs_PVM$fdr) 
  enr.specificPVML_vs_PVM$OR.odds.ratio <- as.numeric(enr.specificPVML_vs_PVM$OR.odds.ratio) 
  enr.specificPVML_vs_PVM$NOM_pval <- as.numeric(enr.specificPVML_vs_PVM$NOM_pval)
  topspecificPVML_vs_PVM <- topPath(fishBP ,"specificPVML_vs_PVM", 0.7, 20)
  enr.specificPVML_vs_PVM <- enr.specificPVML_vs_PVM[match(topspecificPVML_vs_PVM, rownames(enr.specificPVML_vs_PVM)),]
  rownames(enr.specificPVML_vs_PVM) <- gsub("GOBP_", "",rownames(enr.specificPVML_vs_PVM)) 
  enr.specificPVML_vs_PVM=enr.specificPVML_vs_PVM[order(enr.specificPVML_vs_PVM$NOM_pval),]
  enr.specificPVML_vs_PVM$comp = "specificPVML_vs_PVM"
  
  enr.specificPVML_vs_NAWM = fishBP$specificPVML_vs_NAWM
  enr.specificPVML_vs_NAWM$pathway = rownames(enr.specificPVML_vs_NAWM)
  enr.specificPVML_vs_NAWM$fdr <- as.numeric(enr.specificPVML_vs_NAWM$fdr) 
  enr.specificPVML_vs_NAWM$OR.odds.ratio <- as.numeric(enr.specificPVML_vs_NAWM$OR.odds.ratio) 
  enr.specificPVML_vs_NAWM$NOM_pval <- as.numeric(enr.specificPVML_vs_NAWM$NOM_pval)
  enr.specificPVML_vs_NAWM=enr.specificPVML_vs_NAWM[order(enr.specificPVML_vs_NAWM$NOM_pval),]
  topspecificPVML_vs_NAWM = subset(enr.specificPVML_vs_NAWM,NOM_pval<=0.1)
  topspecificPVML_vs_NAWM = topspecificPVML_vs_NAWM$pathway
  enr.specificPVML_vs_NAWM <- enr.specificPVML_vs_NAWM[match(topspecificPVML_vs_NAWM, rownames(enr.specificPVML_vs_NAWM)),]
  rownames(enr.specificPVML_vs_NAWM) <- gsub("GOBP_", "",rownames(enr.specificPVML_vs_NAWM)) 
  enr.specificPVML_vs_NAWM=enr.specificPVML_vs_NAWM[order(enr.specificPVML_vs_NAWM$NOM_pval),]
  enr.specificPVML_vs_NAWM$comp = "specificPVML_vs_NAWM"
  enr.specificPVML_vs_NAWM <- head(enr.specificPVML_vs_NAWM,20)
  
  enr.specificNAGM_vs_GM = fishBP$specificNAGM_vs_GM
  enr.specificNAGM_vs_GM$pathway = rownames(enr.specificNAGM_vs_GM)
  enr.specificNAGM_vs_GM$fdr <- as.numeric(enr.specificNAGM_vs_GM$fdr) 
  enr.specificNAGM_vs_GM$OR.odds.ratio <- as.numeric(enr.specificNAGM_vs_GM$OR.odds.ratio) 
  enr.specificNAGM_vs_GM$NOM_pval <- as.numeric(enr.specificNAGM_vs_GM$NOM_pval)
  enr.specificNAGM_vs_GM=enr.specificNAGM_vs_GM[order(enr.specificNAGM_vs_GM$NOM_pval),]
  topspecificNAGM_vs_GM = subset(enr.specificNAGM_vs_GM,NOM_pval<=0.1)
  topspecificNAGM_vs_GM = topspecificNAGM_vs_GM$pathway
  enr.specificNAGM_vs_GM <- enr.specificNAGM_vs_GM[match(topspecificNAGM_vs_GM, rownames(enr.specificNAGM_vs_GM)),]
  rownames(enr.specificNAGM_vs_GM) <- gsub("GOBP_", "",rownames(enr.specificNAGM_vs_GM)) 
  enr.specificNAGM_vs_GM=enr.specificNAGM_vs_GM[order(enr.specificNAGM_vs_GM$NOM_pval),]
  enr.specificNAGM_vs_GM$comp = "specificNAGM_vs_GM"
  enr.specificNAGM_vs_GM <- head(enr.specificNAGM_vs_GM,20)
  
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = fishBP$commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$pathway = rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM)
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$fdr <- as.numeric(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$fdr) 
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$OR.odds.ratio <- as.numeric(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$OR.odds.ratio) 
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval <- as.numeric(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval)
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM=enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM[order(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval),]
  topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = subset(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM,NOM_pval<=0.1)
  topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$pathway
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM <- enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM[match(topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM, rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM)),]
  rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM) <- gsub("GOBP_", "",rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM)) 
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM=enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM[order(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval),]
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$comp = "common PVML_vs_NAWM & PVML_vs_WML & PVML_vs_PVM"  
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM <- head(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM,20)

Specific WML_vs_CL

ggplot(enr.specificWML_vs_CL, aes(x=reorder(rownames(enr.specificWML_vs_CL),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "WML_vs_CL") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Specific PVML_vs_PVM

ggplot(enr.specificPVML_vs_PVM, aes(x=reorder(rownames(enr.specificPVML_vs_PVM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "PVML_vs_PVM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Specific PVML_vs_NAWM

ggplot(enr.specificPVML_vs_NAWM, aes(x=reorder(rownames(enr.specificPVML_vs_NAWM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "PVML_vs_NAWM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Specific NAGM_vs_GM

ggplot(enr.specificNAGM_vs_GM, aes(x=reorder(rownames(enr.specificNAGM_vs_GM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "NAGM_vs_GM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

common PVML_vs_NAWM & PVML_vs_WML & PVML_vs_PVM

ggplot(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM, aes(x=reorder(rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "common PVML_vs_NAWM & PVML_vs_WML & PVML_vs_PVM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Molecular Function

  enr.specificWML_vs_CL = fishMF$specificWML_vs_CL
  enr.specificWML_vs_CL$pathway = rownames(enr.specificWML_vs_CL)
  enr.specificWML_vs_CL$fdr <- as.numeric(enr.specificWML_vs_CL$fdr) 
  enr.specificWML_vs_CL$OR.odds.ratio <- as.numeric(enr.specificWML_vs_CL$OR.odds.ratio) 
  enr.specificWML_vs_CL$NOM_pval <- as.numeric(enr.specificWML_vs_CL$NOM_pval)
  enr.specificWML_vs_CL=enr.specificWML_vs_CL[order(enr.specificWML_vs_CL$NOM_pval),]
  topspecificWML_vs_CL = subset(enr.specificWML_vs_CL,NOM_pval<=0.1)
  topspecificWML_vs_CL = topspecificWML_vs_CL$pathway
  enr.specificWML_vs_CL <- enr.specificWML_vs_CL[match(topspecificWML_vs_CL, rownames(enr.specificWML_vs_CL)),]
  rownames(enr.specificWML_vs_CL) <- gsub("GOMF_", "",rownames(enr.specificWML_vs_CL)) 
  enr.specificWML_vs_CL=enr.specificWML_vs_CL[order(enr.specificWML_vs_CL$NOM_pval),]
  enr.specificWML_vs_CL$comp = "specificWML_vs_CL"
  enr.specificWML_vs_CL <- head(enr.specificWML_vs_CL,20)
  
  enr.specificPVML_vs_PVM = fishMF$specificPVML_vs_PVM
  enr.specificPVML_vs_PVM$pathway = rownames(enr.specificPVML_vs_PVM)
  enr.specificPVML_vs_PVM$fdr <- as.numeric(enr.specificPVML_vs_PVM$fdr) 
  enr.specificPVML_vs_PVM$OR.odds.ratio <- as.numeric(enr.specificPVML_vs_PVM$OR.odds.ratio) 
  enr.specificPVML_vs_PVM$NOM_pval <- as.numeric(enr.specificPVML_vs_PVM$NOM_pval)
  enr.specificPVML_vs_PVM=enr.specificPVML_vs_PVM[order(enr.specificPVML_vs_PVM$NOM_pval),]
  topspecificPVML_vs_PVM = subset(enr.specificPVML_vs_PVM,NOM_pval<=0.1)
  topspecificPVML_vs_PVM = topspecificPVML_vs_PVM$pathway
  enr.specificPVML_vs_PVM <- enr.specificPVML_vs_PVM[match(topspecificPVML_vs_PVM, rownames(enr.specificPVML_vs_PVM)),]
  rownames(enr.specificPVML_vs_PVM) <- gsub("GOMF_", "",rownames(enr.specificPVML_vs_PVM)) 
  enr.specificPVML_vs_PVM=enr.specificPVML_vs_PVM[order(enr.specificPVML_vs_PVM$NOM_pval),]
  enr.specificPVML_vs_PVM$comp = "specificPVML_vs_PVM"
  enr.specificPVML_vs_PVM <- head(enr.specificPVML_vs_PVM,20)
  
  enr.specificPVML_vs_NAWM = fishMF$specificPVML_vs_NAWM
  enr.specificPVML_vs_NAWM$pathway = rownames(enr.specificPVML_vs_NAWM)
  enr.specificPVML_vs_NAWM$fdr <- as.numeric(enr.specificPVML_vs_NAWM$fdr) 
  enr.specificPVML_vs_NAWM$OR.odds.ratio <- as.numeric(enr.specificPVML_vs_NAWM$OR.odds.ratio) 
  enr.specificPVML_vs_NAWM$NOM_pval <- as.numeric(enr.specificPVML_vs_NAWM$NOM_pval)
  enr.specificPVML_vs_NAWM=enr.specificPVML_vs_NAWM[order(enr.specificPVML_vs_NAWM$NOM_pval),]
  topspecificPVML_vs_NAWM = subset(enr.specificPVML_vs_NAWM,NOM_pval<=0.1)
  topspecificPVML_vs_NAWM = topspecificPVML_vs_NAWM$pathway
  enr.specificPVML_vs_NAWM <- enr.specificPVML_vs_NAWM[match(topspecificPVML_vs_NAWM, rownames(enr.specificPVML_vs_NAWM)),]
  rownames(enr.specificPVML_vs_NAWM) <- gsub("GOMF_", "",rownames(enr.specificPVML_vs_NAWM)) 
  enr.specificPVML_vs_NAWM=enr.specificPVML_vs_NAWM[order(enr.specificPVML_vs_NAWM$NOM_pval),]
  enr.specificPVML_vs_NAWM$comp = "specificPVML_vs_NAWM"
  enr.specificPVML_vs_NAWM <- head(enr.specificPVML_vs_NAWM,20)
  
  enr.specificNAGM_vs_GM = fishMF$specificNAGM_vs_GM
  enr.specificNAGM_vs_GM$pathway = rownames(enr.specificNAGM_vs_GM)
  enr.specificNAGM_vs_GM$fdr <- as.numeric(enr.specificNAGM_vs_GM$fdr) 
  enr.specificNAGM_vs_GM$OR.odds.ratio <- as.numeric(enr.specificNAGM_vs_GM$OR.odds.ratio) 
  enr.specificNAGM_vs_GM$NOM_pval <- as.numeric(enr.specificNAGM_vs_GM$NOM_pval)
  enr.specificNAGM_vs_GM=enr.specificNAGM_vs_GM[order(enr.specificNAGM_vs_GM$NOM_pval),]
  topspecificNAGM_vs_GM = subset(enr.specificNAGM_vs_GM,NOM_pval<=0.1)
  topspecificNAGM_vs_GM = topspecificNAGM_vs_GM$pathway
  enr.specificNAGM_vs_GM <- enr.specificNAGM_vs_GM[match(topspecificNAGM_vs_GM, rownames(enr.specificNAGM_vs_GM)),]
  rownames(enr.specificNAGM_vs_GM) <- gsub("GOMF_", "",rownames(enr.specificNAGM_vs_GM)) 
  enr.specificNAGM_vs_GM=enr.specificNAGM_vs_GM[order(enr.specificNAGM_vs_GM$NOM_pval),]
  enr.specificNAGM_vs_GM$comp = "specificNAGM_vs_GM"
  enr.specificNAGM_vs_GM <- head(enr.specificNAGM_vs_GM,20)

  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = fishMF$commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$pathway = rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM)
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$fdr <- as.numeric(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$fdr) 
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$OR.odds.ratio <- as.numeric(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$OR.odds.ratio) 
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval <- as.numeric(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval)
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM=enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM[order(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval),]
  topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = subset(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM,NOM_pval<=0.1)
  topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM = topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$pathway
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM <- enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM[match(topcommonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM, rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM)),]
  rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM) <- gsub("GOMF_", "",rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM)) 
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM=enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM[order(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$NOM_pval),]
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM$comp = "common PVML_vs_NAWM & PVML_vs_WML & PVML_vs_PVM"
  enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM <- head(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM,20)

Specific WML_vs_CL

ggplot(enr.specificWML_vs_CL, aes(x=reorder(rownames(enr.specificWML_vs_CL),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "WML_vs_CL") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Specific PVML_vs_PVM

ggplot(enr.specificPVML_vs_PVM, aes(x=reorder(rownames(enr.specificPVML_vs_PVM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "PVML_vs_PVM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Specific PVML_vs_NAWM

ggplot(enr.specificPVML_vs_NAWM, aes(x=reorder(rownames(enr.specificPVML_vs_NAWM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "PVML_vs_NAWM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

Specific NAGM_vs_GM

ggplot(enr.specificNAGM_vs_GM, aes(x=reorder(rownames(enr.specificNAGM_vs_GM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "NAGM_vs_GM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))

common PVML_vs_NAWM & PVML_vs_WML & PVML_vs_PVM

ggplot(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM, aes(x=reorder(rownames(enr.commonPVML_vs_NAWM.PVML_vs_WML.PVML_vs_PVM),-rank(NOM_pval)), y=-log10(NOM_pval),fill=-log10(NOM_pval))) + geom_col(width=.65)+ scale_fill_gradient(low="#F6BDC0",high="#DC1C13") +theme(axis.text=element_text(size=8)) + scale_x_discrete(label = function(x) stringr::str_wrap(gsub("_", " ",x), width = 30)) + labs(title = "common PVML_vs_NAWM & PVML_vs_WML & PVML_vs_PVM") + xlab("") + coord_flip() + theme_bw()+ theme(axis.text = element_text(face="bold"))